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Abstract We describe how to apply the recursive Green's 
function method to the computation of electronic transport 
properties of graphene sheets and nanoribbons in the linear 
response regime. This method allows for an amenable in- 
clusion of several disorder mechanisms at the microscopic 
level, as well as inhomogeneous gating, finite temperature, 
and, to some extend, dephasing. We present algorithms for 
computing the conductance, density of states, and current 
densities for armchair and zigzag atomic edge alignments. 
Several numerical results are presented to illustrate the use- 
fulness of the method. 

Keywords electronic transport • recursive Green's function 
method • graphene nanoribbons 
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1 Introduction 

Since graphene was first produced, several synthesis strate- 
gies have been put forward. Significant progress has been 
made to produce better quality samples with the goal of im- 
proving their transport properties. Despite the enormous ef- 
fort, we are still very far from reaching the perfect ballistic 
regime and disorder always plays a central role, particularly 
in electronic transport. Disorder appears in several different 
forms, being either local (such as lattice defects, edge ir- 
regularities, and surface adsorbates) or long ranged (such as 
charge impurities trapped in the substrate or ripples due to 
substrate roughness) |[T]. 
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Several theoretical methods have been developed to de- 
scribe electronic transport in disordered graphene [0- The 
effective low-energy Dirac Hamiltonian, derived from band- 
structure theory, combined with a standard diagrammatic 
perturbative expansion is an excellent analytical tool for giv- 
ing us insight into the properties of disordered graphene B 
m. However, it has {kpi)^^ as a small expansion parameter, 
where kf stands for the Fermi wave number and i is the elec- 
tron mean free path. Thus, it describes well the conductivity 
in graphene at high doping, but becomes of limited use when 
one is interested in the physics close to the charge neutrality 
point, where kf£ ^ 1. Theoretical investigations of the that 
regime require instead the use of numerical methods. 

Most numerical methods employed to study the trans- 
port properties of disordered graphene use an atomistic ba- 
sis IJl. The few exceptions are tailor-made methods to deal 
with long-range disorder, where either a momentum repre- 
sentation ||5]|6l or discretized version of the Dirac equation 
ITUH are used within a single-valley approximation. 

For many applications, one is interested in the two- or 
multiple-probe conductance. For the conductance, differently 
from the conductivity, geometry plays an important role. 
The recursive Green's function (RGF) method ^ became 
the standard tool to compute transport properties in this case. 
The method is very reliable, computationally efficient, and 
allows for a parallel implementation [jCll . It can model arbi- 
trary geometries and efficiently addresses a variety of scat- 
tering processes within the single-particle approximation. 
The goal of this paper is to show how to compute electronic 
transport properties of graphene samples within the tight- 
binding approximation using the RGF method. The key el- 
ement is an efficient algorithm for evaluating the single- 
particle Green's function of sheets or ribbons. 

The recursive method was developed by Thouless and 
Kirkpatrick ||9l for computing the linear electronic conduc- 
tance of Unear atomic chains in the presence of on-site dis- 
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order. The method was later generaHzed to two-dimensional 
systems in the "slice" formulation, which is the form most 
used nowadays ifTTIl . Variations of the method have been in- 
troduced in the Uterature to treat three-dimensional lITZl and 
multi-probe systems lfT3l with arbitrary geometries, see e.g. 
Ref. lfT4l . 

We note that other efficient, atomistic methods have been 
employed in recent years to study electronic transport in 
mesoscopic systems: For instance, the wave-packet time evo- 
lution II15II161 . the kernel polynomial expansion II17II18I . and 
the continued fraction expansion lfT9l . to name a few. Re- 
cently, an alternative method to compute transport of bal- 
listic graphene junctions, particularly effective when strong 
magnetic fields are presented, was introduced [l20l . How- 
ever, for most practitioners, the RGF remains the best method 
for tackling large-scale but finite-size problems where quan- 
tum coherence and disorder are present simultaneously. 

This paper does not attempt to be a comprehensive re- 
view of the recursive method, but rather a self-contained de- 
scription that gives to interested readers, yet unfamiliar with 
the RGF method, all the basic material necessary to imple- 
ment a calculation on their own. For that purpose, we briefly 
present some standard material covered in textbooks iBTI 
I22II . discuss some more advanced issues which are found 
scattered in the literature, and present original developments 
tailor-made for graphene. 

This paper is organized as follows. We begin by quickly 
reviewing some fundamental relations of electronic trans- 
port theory and by providing the essential formulation of the 
method. In Sec.[3]we present the recursive Green's function 
method. The method requires as input the surface Green's 
functions of the electronic leads, taken at the lead-device in- 
terface. In Sec. |4] we describe how to compute the surface 
Green's function of semi-infinite lattices that play the role 
of leads. Next, we present an efficient discretization scheme 
to implement the recursive method for graphene sheets and 
nanoribbons. In Sec.|6]we show how to evaluate quantities 
such as the local density of states and local current densities. 
Very often, one is interested in cases where the coherence 
length £^ is comparable to the system size L. For such situ- 
ations, it is possible to account for dephasing using the phe- 
nomenological voltage probe model, as described in Sec.|7] 
In the context of graphene, the main application of the RGF 
method is the study of disorder effects in electronic trans- 
port. We discuss the main kinds of disorder and show how 
to account for them in Sec. |8] We conclude by presenting 
a number of numerical results that illustrate the method in 
Sec.|9] 



dimensional systems. The linear dc conductance is computed 
using the exact single-particle retarded Green's function that 
connects the source and drain leads, in conjunction with ei- 
ther the Landauer ||231 or the Caroh formula [l24l . 

For a two-probe setup, as illustrated by Fig. [T] the zero- 
temperature linear conductance is given by the Landauer 
formula 



XT'.['''] 



(1) 



Here, t{t') is the transmission matrix across the system from 
left to right (right to left) and r(r') is the reflection matrix at 
the left-hand (right-hand) side. The factor of 2 stands for 
spin degeneracy and the trace is taken over the propagating 
modes at the left and right leads. The transmission matrix 
can be obtained from the S matrix. 



S = 



rt' 
t r' 



(2) 



which is given by II25II26I 

Sah{E) = - dab + iTl^/VaVb 

X J dyci J dypX*{yci)G'{y^,yp;E)Xb{yp), 



(3) 



where v^ and Xc{yp) ^"e, respectively, the longitudinal prop- 
agation velocity and its transverse wave function in the prop- 
agating channel c of lead p (either on the left-hand or right- 
hand side). The integrations run over the contact regions at 
the right and left terminations of the graphene sheet (see Fig. 
[T]). The key element in Eq. (O is G''{yci,yp',E), the retarded 
Green's function corresponding to an electron with energy 
E propagating from positions yp to y^. 
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Fig. 1 Typical two-probe scheme used in the numerical calculations: 
The sample is described by G(yg,yp;E). The perfect leads can be ac- 
counted for by either propagating mode wave functions Xc and their 
density of states Pc (Landauer formula) or by level widths F (Caroli 
formula). 



2 Elements of Linear Mesoscopic Transport 

In this Section we review the key elements necessary to 
implement the recursive Green's function method for two- 



The recursive Green's function method reviewed in this 
paper is an efficient tool to compute the scattering proper- 
ties of noninteracting electrons described by a tight-binding 
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Hamiltonian of the form 

H=-Y^{t.,\i){j\+H.c)+Y^Vi\i){i\ 



(4) 



In graphene, the hopping terms f,j typically connect only 
nearest neighbor sites / and j of a honeycomb lattice (single- 
orbital approximation), although it is straightforward to in- 
clude next-to-nearest hopping terms as well, if required. The 
model can account for an external magnetic field by a suit- 
able modification of the hopping terms, as shown in Ap- 
pendix |B] Here, V, stands for a local potential due to gating 
or disorder 

For simplicity, here we consider only tight-binding mod- 
els with orthogonal orbitals and nearest-neighbor hopping 
terms. The computational cost of the RGF method scales as 
A^ X M^, where A^ is the number of slices and M is the typi- 
cal number of sites in a given slice, see Fig. [1] The RGF im- 
plementation scheme presented in this paper is particularly 
recommended when the system's translational invariance is 
broken by disorder and/or an irregular geometry. For sys- 
tems with translational invariance at the transverse direction, 
it can be advantageous to work in the fc-space and use alter- 
native hybrid RGF implementations (such as in Ref. liZTl ) or 
other methods, e.g. Ref. 1201 . 

In the tight-binding basis, Eq. (O reads 



SalAE) 



Oab - 



ih- 



ao 



■LLx*aii)G',p{iJ)XbiJ), (5) 



'^pj^i 



where the sums run over the sites at the contacts p and q 
where the propagating channels a and b are defined, respec- 
tively. Notice that in two spatial dimensions 



Xaiyp) 



1 



and 



Xa{i) 



1 



[i\Xa) 






(6) 



(7) 



where gq is the lattice constant. Let us consider the case 
where the leads are modeled by semi-infinite square lattices. 
One can then introduce the level widths 11211 

rAi:i')^LXc.ii)—Xa{i'). (8) 

It is straightforward to show that, in this case, 

i.i'eL 



where 

=T = Trj rLGl^FRGgi^ 



(9) 



(10) 



Here, the subscript in the trace indicates whether the sums 
run over channels (c) or sites (s). Depending on the author, 
the expression on the rh.s. of Eq. (|9]i is called either Caroli 
ll24l or Meir-Wingreen ll28l conductance formula. 

This demonstration of the equivalence between the Lan- 
dauer and Caroli formulas relies on the Fisher and Lee 5- 
matrix and on an expression for Fp which is only suitable for 
a square lattice. This derivation is simple and to some extend 
non rigorous but captures the essential elements that will 
be discussed in what follows, namely, the Green's functions 
G^'l and the decay width matrices Frj^. There are several 
ways to show that ^ holds in general in the linear response 
regime, see e.g. Ref. Il29l . 

When the full S matrix is known, it is possible to obtain 
the global density of states through the Wigner time delay 

1, namely. 






(11) 



where the derivative of S with respect to the energy can be 
done numerically. The computation of ( fTTT i is significantly 
less expensive than evaluating p{E) through the standard 
expression, namely. 



p(£) = --Im Tr/GV^) 
n 



(12) 



but it requires the knowledge of the explicit form of the lead 
wave functions Xa{i)- Note that, in Eq. (fT2] i. the trace is 
taken over all sites of the graphene sample. 

The Fano factor is another quantity of interest ll30l . It 
can be evaluated through the expression 



F = 1- 



Tr^ft 



Tr,[f+f] 



(13) 



Notice that one can define left-to-right and right-to-left Fano 
factors, as in the case for the conductance, by switching the 
matrix t with t'. The Fano factor can be computed without an 
explicit knowledge of the wave functions Xa{i) by noticing 
that Eq. (fT3b can be recast as 



1 



Trj FLG'[^ifFRG"i^[^FLGli^FRG"[([_ 



sr 



(14) 



Let us present the same basic expression in a more suit- 
able form for the recursive calculations. We begin by writing 
the Caroli formula for the transmission probability at a given 
energy E using the slice indexing. 



.'J - Tr, [rLGojv+il-E) JrGn+i,o('E')] : 



(15) 



where G'" are the retarded and advanced Green's functions 
across the system (see Fig.|2]for a definition of the subscripts 
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in terms of slice numbers). These Green's functions are ma- 
trices whose rank is defined by the number of sites in the 
slicesljThe level width matrices are given by the expression 



Dyson formulas for an exact Green's function (for a deriva- 
tion of these formulas, see Refs. 12111311 ). 



rL.R = i Ei^{E)-L[';^{E) , 



(16) 



where the retarded surface self-energies of the leads read 

Z[{E)^iiLg'L{E)ul and i:^(£) = 4^^(£)mk. (17) 

Notice that the retarded Green's functions g£ and g^ are de- 
fined at the surface of the left and right leads, respectively, 
when the leads are decoupled from the system. They obey 
the self-consistent equations 

[E + iO+-hL-E[{E)]gl{E)=I and 

[E + iQ+-hR-E'j,{E)]g'^{E)^I, (18) 

where h^ and Iir are the Hamiltonians of isolated, individ- 
ual slices in the left and right leads, respectively. The con- 
nection matrices m/, and ur are defined to run from left-to- 
right and are assumed uniform inside the leads. If the leads 
are identical, then g^ ~ ^£, hi = Ijr, mr — uj^, E^{E) = 
E[{E), and Fr ~ r^. Notice that, in general, the coupling 
matrices are Hermitian, (Ir^l) = J/j,l, while (C) ~ G". 
These two properties guarantee that the transmission prob- 
ability computed with Eq. (fTSl l is always real. Moreover, 
since these coupling matrices are also positive by their def- 
inition in Eq. ( fT6b [notice that the imaginary part of the 
retarded self-energy is negative if we adopt Eq. (118^1. one 
can show that the transmission probability is positive, as it 
should be. 



"l Uo,i Ui,2 
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Fig. 2 Slicing sciieme. The central rectangle containing the dark strips 
(slices) represents the bulk of the sample 



3 Recursive Green's Functions 

We now present the recursive Green's function method, a 
very efficient way to compute the Green's functions that ap- 
pear in Eq. ( fTSb . We begin by introducing two equivalent 

' The number of sites per slice does not need to be equal for all 
slices. 



G = d'>^ + GVG''''\ 



(19) 
(20) 



where G'"' represents the "unperturbed" Green's function 
and V the perturbation. We use these expressions to obtain 
recursive relations for the exact Green's function of a quasi- 
one-dimensional system coupled to leads. The basic idea is 
to break up the system into independent parts (leads and 
slices) and associate to these parts "unperturbed" Green's 
functions G'"' . The hoping matrix elements connecting those 
parts are then selectively built into the perturbation V. By 
choosing the connecting matrix elements and applying Eqs. 
iT% and (l20b iudiciouslv. we can build the full Green's func- 
tion G slice by slice. 

Our presentation is specialized to the case of two-probe 
conductance, see Fig. [1] It is necessary to derive several 
intermediate recurrence formulas before obtaining expres- 
sions for the exact Green's function. We first run the recur- 
rence from left to right, generating a family of Green's func- 
tions G^. Thus, at every step, Eq. (fT9T l is employed using 
a different choice for G^^' and V. We repeat the procedure 
from right to left, generating another set of functions G^. Fi- 
nally, we join these two families to obtain the exact G for the 
whole system. In this way, all parts and connecting matrix 
elements are used (and never double counted). 

The system is broken into A^ thin slices, each one carry- 
ing a maximum of M sites or cells, as show in Fig. |2] The 
slices with numbers lower than 1 or larger than A^ represent 
the left and right leads, respectively. The corresponding re- 
tarded surface Green's functions (when the leads are decou- 
pled from the system) are denoted by giiE) and gR{E), as 
noted earlier These Green's functions are computed sepa- 
rately and before the recurrence procedure (see Sec. 2]. The 
retarded Green's function of the isolated «th slice in the sys- 
tem, g„ (£) = (£■ — h„ + /0+) , does not need to be individ- 
ually evaluated before the recursive calculations. Here, h„ 
denotes the Hamiltonian of the isolated «th slice. 

Neighboring slices within the sample are connected to 
each other through the matrices Un-i,n (left to right) and 
[U„-i,„]'^ = C/,,,,,-1 (right to left), with n=l,...,N. The first 
and last slices in the system are connected to their nearest 
neighboring slices in the leads through the coupling matri- 
ces Uo,i and LV,a?+i. The matrix elements of these matrices 
are the tight-binding hopping ampUtudes connecting sites at 
different slices. 

Here, we assumed that the matrices U only connect nea- 
rest-neighbor slices. For tight-binding models that include 
next-nearest hopping terms, one can still use this algorithm 
by doubling the "width" of the unit slices, which slows down 
the computation by a factor 2^. It also is possible to deal with 
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next-nearest hopping terms and incur in a smaller slowdown 
factor by properly modifying the recursive method ll32l . 

We use subscripts to denote longitudinal spatial indices 
(except for gi, gR, and g„). Thus, Gn.m{E) is the matrix 
Green's function connecting the n and m slices. Sites in- 
dices are shown as a pair of variables: G„,m{j,j') denotes 
the Green's function connecting site j in the nth slice to site 
/ in the mth slice. Hereafter, we will drop the energy vari- 
able E (since scattering is assumed elastic, E is conserved 
throughout the system). 



3.1 Connection to leads 

For the two-terminal setup we address here, the sample (cen- 
tral region) is coupled to a left lead L and to a right lead R. 
In the following we show how to built the Green's function 
that describes this coupling. 

We begin finding the Green's function G^. We recall that 
the rightmost slice of the left lead if denoted by 0. Our goal 
is to obtain Gq „ and G^ „ in order to describe electron propa- 
gation in the sample when the left lead is taken into account. 
The reason will become clear when we reach Sec. 13.41 

The first step is to incorporate the « = 1 slice to the left 
contact Green's function g^. This kind of operation is re- 
peated throughout the method and therefore we present it 
in detail. For this purpose, we introduce the kets |0) and |1) 
which represent the states where electrons are found in slices 
n = and « = 1, respectively. The "unperturbed" Green's 
function in this case is G*"' = |0)gL(0| + | l)gi (1 1 while V = 
|0)f/(),i(l| + |l)'^i,o(0| is the perturbation that connects the 
n ~ I slice to the left lead. Then, using Eq. ( fT9] l. we obtain 

(1|G^|1) = (IIgC^II) + ^ (1|g"'V>('w|^I'w'>('w'|G^|1> 



(llG^^'ll) + (l|G"^'|l)(l|y|0)(0|G^|l) (21) 



and 



(0|G^|1) = (0|G("'|1)+ ^ {Q\G°\m){m\V\m'){m\G^\l) 



70|r\\ /r\lT/l 1 \ /I l/^^l 



= (0|G"|0)(0|y|l)(l|G^|l). (22) 

Adopting the more compact notation {n\G^\m) = G^^, we 
drop the bras and kets and can rewrite these equation as 



Gu=^i+^if/i,oGo,i 
and 

Therefore, 

Gfj ^{I-giUi.QgLUo.iy^ gi. 

Now, since gi = {E — hi)^ ,we can write 

Gf, =(£-/^l-t/l,ogLf/o,I)"'• 



(23) 



(24) 



(25) 



(26) 



Notice that this Green's function takes into account the cou- 
pling of the first slice with the left lead, but has no informa- 
tion about the rest of the system or the right lead. 

It is important to remark that we neglected the infinitesi- 
mal imaginary part in ^i because we assumed that the "self- 
energy" term in Eq. (l26l l brings its ownfinite imaginary part. 

We proceed analogously in order to connect the last slice 
to the right lead. Choosing G^"^ =gR+gN and V ~ Un.n+i, 
we have 



Gn.N — gN + gN Un.N+ 1 G^+ 1 ^N 

and 



Gn+ l,N — gR ^N+ 1 .N G^jj . 



(27) 



(28) 



(Note that the slice indices for the right Green's functions 
run opposite to those in the left Green's functions.) There- 
fore, 



Gn.N — (^ ^gN Un.N+ 1 gR Un+ \.n) gN- 

Again, since gN = [E — h^y , we can write 



Gn.n — {E — hN - UN.N+\gR Un+\.n) 



-I 



(29) 



(30) 



The Green's function G\ j (or Gf/^) describes all single- 
electron processes that begin and end that on the « = 1 (or 
n = A^) slice, taking into account all possible number of in- 
cursions in and out of the left (or right) lead. It does not yet 
take into account incursions into the bulk of the system. 



3.2 Left Green's functions 

With Gf J in hand, we can evaluate the next successive A^ — 1 
left Green's functions by using a recurrence formula analo- 
gous to Eq. ( |26] |. To derive such formula, we choose G'"' = 
Gn-i.n-i ^"d ^ = U„-i^„ + f/„,„-i. Applying Eq. (O, we 
write 



-1 



G„.„ — U -§«^n,n-l G„_i,„_i [/„-!,„) gn, 



-1 



with n = 2,...,N. Using g„ = (£ — h„) , we obtain 

G„.„ = [E — h„ ~ U„_„^\G„_i„_iUn-\,n) 



(31) 



(32) 



This formula is accompanied by another one, which con- 
nects the left-most slice (the surface slice of the left lead) 
with the nth one. 



Gq.„ — Go,„_l t/n-KnG„^„. 



(33) 



Note that A^ inversions are necessary to arrive at the A^th 
slice. Each inversion requires 0{M^) operations. Thus, the 
complexity of the calculation scales as N x M^. 



Caio H. Lewenkopf, Eduardo R. Mucciolo 



3.3 Right Green's functions 

Similarly to left case, for the right Green's functions, using 
Eq. ( l30l l and starting from the Mh slice, we find that 

G„,„= {l-gnU„,„+lG„^l,j^lUn+l,n) gn, (34) 

with« =A^— 1,...,1. Substituting §„ = {E — h„y , we ob- 
tain 

^n.n = [^ ^ ^n — UnM+\ G„+i^„_|_i t/„+l,«) • (35) 

Also, 

^N+Ln = G^W+l.K+l^n+l.nGnn. (36) 

Again, A^ additional inversions have to be performed in or- 
der to arrive at slice the first slice (n = 1), with an overall 
computation cost 0(N x M^). 



3.4 Full Green's functions 

Suppose one arrives at the n slice by either a left or right 
sweep (1 < n < A^). To obtain the exact full Green's function 
of the system we use again Eq. ( fT9l ) assuming G'"^ = g„ + 

G^Ln-l + Gn+ln+1' with V = Un-l,n + U„,„-l + f/„,„+l + 

U„+i.n- As a result, we find 

Gn,n = gn+gn {Un.,n-\ G„-\^n + t^H,n+l G„+i^„) , (37) 



and 
Thus, 



Gn,n — t ^ gn {U„,„-l G„_i „_; C/„_l,„ 

+U„.„+l G„_|_i „+i Un+l,„) 

and since g„ = {E — h„)^ , we obtain 



(38) 



(39) 



-1 



gn, (40) 



Un,n+\ G„+i „+i Un+\,n 

together with 

Go,n = GQi,_iUn^i^„Gn,n 

and 

Gn+Im = G^^ijj^i U„^i„G„_n- 



(41) 



(42) 



(43) 



Note that in order to compute G„.„ and Gat+i,,,, we need 
to keep track of G^„ and G^„ [obtained recursively from 



Eqs. ( l32b and ( |35] ), respectively], as well as Gq,, and G^_|_j „ 
[which follow from Eqs. ( |33] | and ( |36] |. respectively]. In or- 
der to obtain G„_i.„ and G„.„+i, we can apply Dyson's equa- 
tion again to a situation where only the «th slice is decou- 
pled, yielding 



Gn.n+l — G„,„ty„„-|-1 G,j_|_[^ 



+ l,n+li 



while 



G„_ 



i,« 



■'n—l.n—l 



^n— \.n ^n.n- 



(44) 



(45) 



These equations are useful for computing the local current 
distribution (Sec.|6]l. 

We note that when computing the exact Green's in Eqs. 
(ITTl i. (l42l i. (031), (l44l i. and ( |45] | we have selectively used each 
matrix C/„ „/ only once. Similarly, at each step, an isolated 
slice Hamiltonian /i,, was used and never repeated. Thus, 
at the end of the calculation of the full Green's function, 
all hoping amplitudes and local potentials of the underlying 
tight-binding model have been used and only once. 

An alternative way to compute full Green's functions, 
which is quite useful if only transmission and reflection ma- 
trices are required, is to close the left (or right) sweep with 
a connection to the right (left) lead: 



1 . For the left sweep, we use Eq. ( l42l i to write 
which is complemented by 



m+i.N+i 



' Un+ 1 .N Gf^fj Unm+ 1 ) 



(46) 



(47) 



obtained from Eq. (HTl i. 
2. For the right sweep, we use instead Eqs. (l47t and (l42t to 
obtain 



and 

Go,o - 






UifiGofi 



[8l 



— Uo,i Gj 1 U[fl 



-1 



(48) 



(49) 



respectively. 

As we will see below, Eqs. ( l46b and ( |49l ) and can be 
used to compute the left-to-right transmission and left re- 
flection matrices, respectively, while Eqs. (|48) and (l47t yield 
the right-to-left transmission and the right reflection matri- 
ces. For systems with inversion symmetry, we expect Goo = 
Gn+i,n+i and Go,n+i = Ga^+io and therefore only one sweep 
(left or right) is necessary for the evaluation of the whole 
scattering matrix. 

For symmetric leads and in the absence of an external 
magnetic field (i.e., time-reversal symmetric systems). 



H. 



,N+l\ 



Gm 



N+lfi 



(50) 
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and only one sweep is necessary. When such conditions are 
not met, one needs both sweeps, namely, from left-to-right 
and from right-to-left, in order to assemble the scattering 
matrix. Moreover, any local observable (such as the local 
density of states or the local current flux), requires Gq^n+i 
as well as G„^„ for all n = 1,...,N. 



3.5 Input Green's functions 

The recurrence relations shown above rely on some input in- 
formation. One needs to define the Green's functions of the 
leads {gi and gR), the Hamiltonian of the isolated slices {h„, 
n = 1 , . . . , A^), and the hopping between slices (the U matri- 
ces) before starting the calculation of the sample's Green's 
function. 

Since both gi and gR are the input Green's functions, it 
is crucial that they have finite imaginary parts. These will 
be dominant and, in practice, we can basically neglect the 
imaginary part when considering g„ (even if E happens to 
coincide with an eigenvalue of an isolated slice, the imag- 
inary parts brought in by coupling to the leads makes the 
Green's function convergent). We will see next how to ob- 
tain contact Green's functions for leads modeled as semi- 
infinite lattices. 



4 Lead Green's Functions 

To satisfactorily model the leads, there are two main physi- 
cal considerations to keep in mind: (a) the source and drain 
leads in typical graphene transport experiments are metallic 
and thus have a high density of states; (b) graphene-metal in- 
terfaces tend to form ohmic contacts. Thus, in the numerical 
simulations, one needs to eliminate or minimize the contact 
resistance associated to band structure mismatch at the con- 
tacts. To address (a), the chemical potential in the leads is 
customarily adjusted to maximize the density of states. To 
address (b), an appropriate lead lattice model, compatible 
with the sample lattice, is chosen to minimize back reflec- 
tions. The leads are usually modeled either by square or hon- 
eycomb semi-infinite lattices |f33| . In certain cases, a com- 
bination of square lattice contacts coupled to semi-infinite 
linear chains are shown to be advantageous to minimize the 
contact resistance |[34l . 

The methods to compute the lead Green's functions can 
be divided into two categories, namely, the iterative recur- 
sive methods ll35l and the eigenchannel decomposition or 
mode matching ones Illlil36l[37l[38l . In the latter, the lead 
Green's function are built with the eigenchannels of the infi- 
nite (translation invariant) corresponding lattice. Except for 
few cases, such as the semi-infinite square lattice discussed 
below, the eigenmodes depend on the longitudinal wave num- 



ber k and the giR cannot be written in closed analytical 
form. 

In this Section we review the eigenchannel decompo- 
sition of giR for semi-infinite square lattices and discuss 
how to couple giR to a graphene device. Next, we present 
the decimation method for semi-infinite lattices 1.351 . De- 
spite the claim that, in general, iterative methods are inferior 
in performance and accuracy than the eigendecomposition 
ones ll36l . the decimation method has the attractive proper- 
ties of being very robust and straightforward, allowing for a 
very amenable implementation. 



4.1 Square lattice leads - analytical approach 

Let us model the contacts by semi-infinite tight-binding squa- 
re lattices ll33l . Let us also set the Fermi energy to £ = 1391 
[33JI . We can shift the energy band of the electronic states in 
the leads by varying a gate potential Viead in the leads 1401 
[39ll33ll . At zero bias and in the absence of inelastic scatter- 
ing, this sets the energy of the electrons propagating through 
the graphene sheet. 

The electron wave function in the semi-infinite square 
lattice is extended along the x-direction and is quantized in 
the transverse direction. Let us consider hard-wall boundary 
conditions at the edges of the strip, namely, j = and j = 
M + 1 . The transverse wave functions are given by 



Zv(y) = 



M+l 



■sin 



M+l 



(51) 



where v = I,- ,M. Associated to each transverse mode 
there are two extended Bloch waves with longitudinal wave 
numbers ±^„, which are real for propagating modes and 
complex for evanescent modes. 

For convenience, we assume l22l a hard-wall boundary 
condition at the left end of the strip (where it connects to the 
right Green's function of the strip), such that (j)^{N—l) ~Q: 



0m(«) 



sin[^(«-jV+l)], 



(52) 



where /i is a longitudinal quantum number When contrasted 
with a continuum model, we can identify jx = kxQx, where 
kx is the longitudinal wave vector By choosing ^^{N — 1) = 
we are not capable of treating the contribution of evanes- 
cent modes. This is usually not a problem in practice, unless 
measurements are done very close to the contacts (i.e., for 
very short systems). 

The dispersion relation for this model is 



TtV 



£'vM = ^lead-2fvCOS^-2fj.COS 

^ \M+1^ 

The velocities of the propagating modes are 
_ oq / dEvfi \ _ 2flofv 



Ti \ dpi 



h 



sin/i. 



(53) 



(54) 
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By setting iivfi = 0, one writes sin/x = ■\/l — cos^/x with 

COS;U = VleadA.r- (fvAA-)cOS [7rv/(M+ 1)]. 

Let us now construct the leads Green's function as in 
Ref. II22I . The general expression 



g{n,j-n',j';E) = dfij^- 
Jo ^fi 



[0m(«)Zv(;)]10m(«')Zv(/)] 



E — Evfi + iO 



(55) 



can be simplified since we are only interest in n = n' = N 
(surface) and £ = 0. It reads 



2 ^ „ f" . sinV 



8LUJ')^-LXviJrXvU') dli 



v=l 



p + qcosil 



(56) 



where 



p = -yiead + 2fvCOS 



TtV 



By integrating over jj., one writes 

M _ 

Gr'(i,/) = Ez;(;)G^'="'(v)zv(/), 

where 



/0+ and q = 2t^. (57) 



(58) 



v=l 






i-Wi- 



(59) 



Equation dSSl l defines the unitary transformation that con- 
verts the Green's function G'*''™'(v) in the channel represen- 
tation into G^™(7,/) in the site representation. 

In Ref. [331, Schomerus argues that these expressions 
can be related to the parameter jj. for graphene semi-infinite 
lattices in the case of armchair edge orientation. He shows 
that G'™'''(v) = -M"''™''''''(Viead)/f?- This correspondence 
works well for armchair orientations because in that case 
propagating modes do not mix; it does not work for zigzag 
or other edge orientations. In Ref. Il39l the authors specu- 
late that transport properties in the presence of bulk disor- 
der should not depend on the graphene orientation, which is 
numerically confirmed in Ref. BUI . Thus, in large-scale nu- 
merical simulations involving bulk disorder, it is worth tak- 
ing advantage of the matching between square-lattice leads 
and graphene armchair leads. 



4.2 Eigenmode decomposition method - square lattice 

The following alternative approach to obtain the surface Gre- 
en's function of a square-lattice lead is helpful since the 
same steps can be repeated for any other lattice and they 
form the basis of the eigendecomposition methods. 



Since adding another slice to a semi-infinite lead should 
not alter its Green's function, we can write that, in the ab- 
sence of any disorder or inhomogeneity. 



Gf,l=^L 



(60) 



when C/o.i = —txt, where t^ > is the horizontal hopping 
matrix element. Then, using Eq. ( l26l l. we obtain the self- 
consistency condition 



gL--(E-lli~t^gL) V 



(61) 



Here, we assumed Uifi = f^/, appropriate for square lattices, 
where / is the identity matrix. In order to solve Eq. (l6ll for 
gL, we notice that since hi is Hermitian, it must be diag- 
onalizable by a unitary transformation: T^/ii T = diag(ev), 
where diag(£v) is a diagonal matrix containing the eigenval- 
ues of /j 1 . Then, 



|L(v) = (£-ev-f2|L(v)) , 



(62) 



where diag(fL(v)) = T^^^r is the diagonal matrix contain- 
ing the eigenvalues of gi. Solving Eq. ( l62b . we find 



8l(v) 



(E- 



2f? 



■± 



'(£-ev)2 1 



4f4 



(63) 



Notice that for jis — £v | < 2tx, this eigenvalue acquires a fi- 
nite imaginary part. Since we are mainly interested in re- 
tarded Green's functions, we choose the negative sign and 
rewrite the equation as 



8liy)- 



(E-Ey) 

(E-£v) 
It} 



l-Jl- 



(£-£v)^ 



/O^ 



1 {E-£yf 



4r? 



|£-ev|>2f,, 

|£-ev|<2f^. 
(64) 



We need now to determine T and {Cv}. For a square 
lattice, this is trivial: h\ describes a one-dimensional chain 
with M sites and vertical hopping matrix elements ty > 
(hard boundary conditions assumed). Then, 



■V7 



2 ■ f ^vj 
sm ' 



-XvU), 



M+l \M+l^ 
where j = l,...,M and v = 1 , . . . ,M, and 

, TtV 
£v = Vlead — 2f,. COS 



M+l 



resulting in 



M 



8lUJ') = Y. Xv{j)8Liv)Xvif)- 



(65) 



(66) 



(67) 



v=l 



In Appendix |Al we show that Eq. ( IFTI ) gives the expected 
steps in the linear conductance. 
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Unfortunately, the above expressions do not directly ap- 
ply to either zigzag or armchair leads, basically because Uq, i 
is not proportional to the identity in these cases and it does 
not commute with the slice Hamiltonian. The solution for 
the general case is nicely presented in Ref. 



4.3 Leads Green's function - decimation method 

We can use the decimation method of Ref. Il35l to evaluate 
numerically the Green's function of leads with arbitrary (but 
translation invariant) lattice structures. 




12 3 " 

Fig. 3 Chain stmcture used in the decimation method. 

The method works as follows. Suppose we want to solve 
(E — H) G{E) —I for the operator G{E), where the Hamil- 
tonian operator// is defined over a chain where each site has 
an arbitrary basis of dimension J but only nearest-neighbor 
inter-site connections exist (see Fig. O. We assume that all 
these connections are represented by operators u (left-to- 
right) and M^ (right- to-left). The local Hamiltonian is iden- 
tical for all sites and denoted by h. Since, by definition, the 
Green's function obeys 



[(£-//)G(/i)]„,„ = 5„,,„, 
it follows that 

{E-h)Go,o =/ + mGi,o, 
{E-h)G\fi = uGi.a + u Go,o, 



{E-h)G„fi = II G„ 



-i.o- 



« G„_ 



1.0, 



(68) 

(69) 
(70) 

(71) 



with « > 1 . Using Eq. ( iTOt . we find that 

G\fi = {E-hy^ [uG2fi + u^Ga,a) , 

which can be combined with Eq. ( |69] | to yield 

{E - h) Go,o = / + M (£ - /j)"' {uGifi + m'^'Go.o) , 

which can be rewritten as 

[E-h-ii{E-h)-'^u''] Go.o=I+u{E-hy^uG2.Q. (74) 

Notice that we can relate G2,o to Gq^o without involving G^q. 



(72) 



(73) 



The same trick can be employed for any value of «. From 
Eq. (fTTT i. we can write 

G„+i.o = {E-hy^ {uG„+2.o + u'G„fl) , (75) 

and 

G„_Lo = {E-hy^ {u'Gn-2.o + uGn.o) ■ (76) 

These equations can be combined with Eq. (fTll to yield 

{E-h)G„fi =Li{E -hy^ [iiG„+2.o + u^G„,o) + 

u'iE-hy^ {uG„fi + u'G„-2fi) , (77) 

which can be rewritten as 
[E-h- II {E - hyhi'^ - u\e - hy^u] G„,o = 

u{E — hyh{G„+2.o + u\E — hy^u^Gn-2.0- (78) 

Equations (l74l i and (ItST i generate a new recursion series 
involving only even sites; 



(£-enGo^o = / + aiG2,o 
(£-ei)G2,o = «iG4,o + /3iGo,o 

{E - £]_) G„fi = ai G„+2fi + Pi Gn-2fi 

with 

ai — u{E — hy u 

j3i = uHE-hy\'^ 

el = h + u{E-hy^u^ 
ei = £\ + u\E-hy^u. 



(79) 
(80) 

(81) 

(82) 
(83) 
(84) 
(85) 



Even though Eqs. ( |79] | to dSTT l involve only even sites (i.e., 
multiples of 2^), they are identical in form to Eqs. (iTOl i and 
dTTT i. Therefore, we can repeat this procedure k times until 
the recursion relations involve only sites that are multiple of 
2^^, namely. 



(£-e;)Go,o=/ + aiG2,o 
(E -£ii)G2kQ = (^kG2k.2.o + hGQ.Q 



iE-£k)G2k.„x) - C'kG2k.[„+x),o- 
with n > 1 and 



-PkG2k.(n_l),Q 



l5k = l5k-iiE-£k-iy'l5k-u 



(86) 
(87) 

(88) 

(89) 
(90) 
(91) 
(92) 



e^' = Ek-i + ak-\{E -£k-\) Pk-u 

Ek = el + Pk-i{E - Ek-\y^ cik-i- 

The decimation can stop when \\ak\\ and \\Pk\\ are suffi- 
ciently small, in which case we can approximate 

Go,o«(£-eir'- (93) 
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and 



Fig. 4 Structure of the square lattice lead. The slice with n = con'e- 
sponds to the "surface" which will be attached to the system. Slices for 
decimation are denoted by the boxes with dashed lir ~ " 



This provides the Green's function for the 
which can then be related to the lead Greer 
andC^^i- 

Note that one needs to add a small po; 
part to E, namely, E ^ E + ir\,m order to g( 
Green's functions. On the practical side, a 
helps to speed up the convergence but spo 
of Gg at the order of r\/t. Provided that r 
an imaginary part to E has Uttle effect on t 
of graphene transport properties for \E\/r\ 
a small energy interval around the charge i 
where £ = 0. 

4.3.1 Square lattice lead 

From Fig.m we see that each decimation sit( 
a regular open vertical chain with J lattice s 

KjJ')=Y.<^ii{j)h{f)Ei., 



and 



-f^jj- 



(95) 



where t is the nearest-neighbor hopping amplitude, j,/ 
\,...,J,ll = l,...,J, 



el{jj') = h{jj') + 2aiijj'). 



(99) 



Figure |5] shows a comparison between the analytical and 
decimation results for the G''{E) of a square lattice of width 
J = M = 6, projected onto the eigenchannel basis. Notice 
the excellent agreement. 
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Fig. 5 Real (a) and imaginary part (b) of the Green function of a semi- 
infinite square lattice in the channel (diagonal) representation, g{n), as 
a function of energy. 7 = 4. The solid lines represent the exact analytic 
result, whereas the squares were numerically obtained through the dec- 
imation method for T] = 10^' and k= 10. 



^MO-)=\/y^sin(^^|, and 



En = —2t cos 



/+! 



From Eqs. (|94] i and ( |95] l, we derive the relations 



£i(;,/) = M;,/) + «i (;,/), 



(96) 



(97) 



(98) 



4.3.2 Honeycomb lattice lead - armchair edges 

For the honeycomb lattice with armchair edges, the elemen- 
tary slice contains / stacked hexagons. Thus, M = 2/+ 1, 
where M is the vertical number of atoms. The slice has the 
structure of a vertical ladder chain (2M atoms). In this case 
there is no simple formula for the eigenstates of h. The ma- 
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(b) 



2M-I ' 


2M 


2M-3 i • 


2M-2 


7 • ( 


S 


5h ■ 


6 


si 1 


4 


; A ' 





Fig. 6 (a) Structure of the honeycomb armchair lead. Note that the 
number of interconnect channels runs from 1 to /, while the number of 
vertical coordinate points in the slice is M = 27 + 1 . The slices for dec- 
imation are defined by the boxed regions, (b) Internal index stmcture 
of the elementary slice. 



(a) 



J I 

M- 



(b) 




Fig. 7 (a) Stmcture of the honeycomb zigzag lead. Note that the num- 
ber of interconnect channels runs from 1 to 7 = M, where M is the 
number of vertical coordinate points in the graphene. The slices for 
decimation are defined by the boxed regions, (b) Internal index struc- 
ture of the elementary slice. 



trix h, which is 2M x 2M dimensional, reads 



-t -t 

~t -t -f 

-f -f -f 

-f -f -f 

-t -t 

-t -t 



V 



/ 



Let us call {^^(m)} its eigenvectors and {E^} its eigenval- 
ues, with jj. = 1 , . . . , 2M. Then, 



2Af 



M;,/)=I'^M(4y)'^M(4/-i)£M- 



(101) 



4.3.3 Honeycomb lattice lead - zigzag edges 

For the zigzag graphene lead, there is a direct correspon- 
dence between slice channels and the vertical indices of atomic 
positions: M ~ J. The elementary slice is a open vertical 
chain of length 2M. Its eigenfunctions and eigenvalues are 



Mn = 



2 . f nul 
sin 



2M+1 \2M+l 

Ttjl 



and 



En = — 2f COS , , , 



(102) 



with j,f — 1 , . . . , 7. Everything else is similar to Sec. 14.3.11 



with / = 1 , . . . ,M and jj. = I,... , 2M. The slice matrix reads 

2M 

KJJ') = L <^M(2y- 1)'/'a:(2/)£m- (103) 

with j,j' — 1, . . . ,M. All other aspects are identical to Sec. 

5 Device Green's Function 

The device Green's function Gij can be used to describe 
nearly all the physics of transport and its calculation is where 
most of the computational time is spent. This is where any 
optimization of the computational method is most welcome, 
particularly in the study of disorder effects in the electronic 
transport when extensive disorder averaging is required. With 
the RGF method one can compute Gjj for a large variety 
of settings. For instance, the graphene sheet can have an 
arbitrary number of layers and different edge orientations , 
namely, zigzag, armchair, and chiral. Also, the tight-binding 
model can include next-nearest neighbor hopping terms in 
addition to the nearest neighbor ones. These elements have 
to be taken into account when choosing the slice unit cell 
employed by the recursive method. The description of dis- 
order modeling is postponed to section |8] 

In this Section we present efficient slicing schemes for 
graphene monolayers with armchair and zigzag edges. 



( 1 00) 5 . 1 Shcing armchair lattices 



For rectangular geometries, the bottleneck of the recursive 
method is the matrix inversion required for adding a new 
slice [see Eq. (|32])1. Thus, it is always important to try to 
minimize the number of sites in the slice. Having this in 
mind, there is a way to mount the armchair slices which 
reduces the number of sites per slice without introducing 
next-to-nearest neighbor hopping. It is based on the lattice 
deformation shown in Fig.|8] 

The relations between lattice size and graphene sheet di- 
mensions for the efficient armchair slicing are 

L V3 fN \ , V3 . ^ ,^ 1 
— = — - --1 +-— and —=M-l. 
ao 2 \2 J o flo 



(104) 
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Fig. 8 Efficient way to slice an armchair graphene ribbon (the so- 
called "pine tree" configuration). Atoms from different sublattices are 
indicated by empty (A) and full (B) circles, uq is the lattice constant. 
The dashed lines indicate vertical slices. 



5.2 Slicing zigzag lattices 

The zigzag geometry is shown in Fig.|9] The real stracture is 
shown on the left-hand side (honeycomb lattice). An equiva- 
lent square lattice with missing vertical bonds is also shown 
(the so-called "brick wall" configuration). For zigzag edge 
ribbons there is no efficient, alternative slicing that mini- 
mizes the number of sites per slice without creating a next- 
to-nearest neighbor connectivity. 

By convention, we assume that the site at the left bot- 
tom comer of the lattice is of type A. In that way, slices with 
even and odd number of sites will alternate as we move hor- 
izontally. Two cases will need to be considered separately: 
M odd and M even. However, both cases share a common 
trend, namely, the presence of dimers and isolated sites (at 
the bottom and/or at the top of the slice). Therefore, find- 
ing the Green's function of an isolated slice is a very simple 
exercise. 

The relations between lattice size and graphene sheet di- 
mensions for the zigzag slicing are 



— = ^r{M-l) + —- and — 
flo 2 6 flo 



A^-1 



y 




z 



(P,+; <?;■+/ 1 <P.;+; 




Fig. 9 Graphene strip with zigzag edges. Atoms from different sub- 
lattices are indicated by empty (A) and full (B) circles, c/q is the lat- 
tice constant. The dashed lines indicate vertical slices and the dotted 
highlights adimer. By convention, we set as (1,1) the coordinates of an 
atom of type A placed at the left bottom comer. The Peierls phases (see 
appendix^ of the "horizontal" hopping matrix elements are indicated. 



6 Evaluating Local Quantities 

In addition to the transmission ,'7, the RGF method can be 
used to calculate other quantities such as the local density of 
states (e.g. BTl ) and the local current density (e.g. 042114311 ). 



6. 1 Local density of states 

The local density of states (LDOS) can be easily evaluated 
if one knows the exact retarded Green's function at a given 
site: 



p{nJ\E) 



Um[Gl„{jJ-E)\. 

K 



(106) 



(105) 



The exact local Green's function is evaluated using Eq. (l4Tb . 
Conductance calculations require just a single sweep through 
the lattice, since the transmission formula (flSt only needs 
G^ evaluated at the slices corresponding to the contacts. In 
contrast, the LDOS demands the calculation of C^ at all lat- 
tice slices of interest, which further increases the compu- 
tational cost linearly with A^, resulting in 0{N^M^). Let us 
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provide an example of a situation where the LDOS plays a 
central role and needs to be evaluated. 

Numerous studies have addressed the possibility of lo- 
cal magnetic moment formation in graphene either due to 
zigzag terminations at the edges in graphene nanoribbons 
or due to vacancies in graphene sheets (for a review, see 
Ref. iSH). The effect can be understood through the Stoner 
mechanism of magnetism, which requires an enhanced LDOS 
when strong electron-electron interaction as present, as orig- 
inally proposed in Ref. B31 using the Density Functional 
Theory (DFT). The tight-binding Hamiltonian, Eq. (|4|i, can 
be modified to reproduce the DFT results by adding a Hub- 
bard mean-field term P31 . namely. 



(107) 



^ = - E ('•AaCj.a + ii-' 



Kj.a 



^y, + [/(«,-,_,) cj„c. 



where the operators c\^ and C; ^ create and annihilate an 
electron of spin projection o at at the site /. As standard, 
(«, (j) is the occupation number and U is the on-site electron- 
electron interaction strength. In graphene, U is usually fitted 
to reproduce the DFT band structure calculations for trans- 
lation invariant systems (for which V,- is constant). As we 
discuss in Sec. [8] by a suitable choice of the parameters V, 
and tij, the Hamiltonian ( I107l i becomes an excellent frame- 
work to model disorder as well. 

In general, the occupation numbers («,,(j) that appear 
in Eq. (I107l i can be obtained from nonequilibrium Green's 
functions ||3TI . In the linear response regime, a simplifica- 
tion allows one to express {ni^a) in terms of an equilibrium 
Green's function, namely. 



(«;,o 



dElm[Gl„.,{jJ-E)]f{E-^). (108) 



where i= {j,n). Note that Eqs. (I107l i and (llOSI l have to be 
solved self-consistently. 

The large number of poles of G''{E) makes impractical 
the integration of the r.h.s. of Eq. (1108b on the real axis, a 
difficulty shared with transport studies on molecular elec- 
tronics, see e.g. Ref. Il46l . An optimized strategy to imple- 
ment an efficient integration along a complex plane contour 
is presented in Refs. 0411471 . 

Another frequent application of LDOS occurs in the eval- 
uation of the local charge density. The latter can be calcu- 
lated from the LDOS through the expression 



nc{nJ;EF) = - 



1 r^p 



m(«j) 



dEp{nJ;E), 



(109) 



where A is the sheet area and jj. {n,j) denotes the local chem- 
ical potential. Note that in «-type regions, Ef > fi and there- 
fore the integral is over positive energies ("electrons"), while 



in /5-type regions, Ef < jj, and the integral is over negative 
energies ("holes"). The local chemical potential is evaluated 
with respect to Ef — 0, 



ll{n,j)^V{nJ), 



(110) 



where the background potential V{n,j) includes the gate 
voltage. 



6.2 Local current density 

The local current density is involves more complex calcula- 
tion than the LDOS. Several methods were developed in the 
literature (see e.g. Ref. B2l ). The bond current between two 
neighboring sites of lattice coordinates {n,j) and {n',f) is 
obtained using the equations-of-motion method for nonequi- 
librium Green's functions (see e.g. ||3T1 ) and reads ll24l 



2e 



'{n.j)^{n'.f) 



dE 



Un,Ajj')G<,^„{f,j) 
~Un'M-.J)G^,n'{hj') 



(111) 



where G^ is the lesser Green's function of the system. |3 
For the tight-binding model with nearest-neighbor hopping, 
there are two situations to consider First, n = «', in which 
case / = y ± 1 and the current is intra-slice. Second, when 
«' = 71 ± 1, the current is inter-slice and / ~ j or / = j ± 1, 
at most. Notice that in the absence of magnetic fields, time- 
reversal symmetry requires U„„'{j,j') = t/„/ „{/ J)- 1" ad- 
dition, in equilibrium conditions, the lesser Green's function 
is a symmetric matrix, leading to a zero bond current, as 
expected. Thus, local currents can only appear through the 
application of a magnetic field (which breaks time-reversal 
symmetry) or when a finite bias voltage between contacts 
exists. 



Equation (llllb requires the calculation of the exact lesser 
Green's function for sites in the bulk of the system. This can 
be done recursively, starting from the equilibrium Green's 
function of the leads (see Sec. 16. 3b . An alternative approach, 
suitable for transport in the Unear regime, is to use the elimi- 
nation technique developed in Ref. B3l . where only retarded 
and advanced Green's functions are required and the en- 
ergy integration is avoided. The computational cost is this 
approach also scales as 0{N^M^). 

For plotting current fields, it is useful to define the local, 
on-site, outgoing vector current as 



I 



"J 



= > a,, 



>(«',/) 



(112) 



^ There are several good textbooks, such as Refs. 131148 1491 , that 
discuss nonequilibrium Green's functions. In particular, Ref. 1311 con- 
cisely covers all required background material. We refer the reader to 
these books for the derivation of expressions involving G^ and related 
functions and further insight into the subject. 
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where the sum is over sites («' , / ) that are nearest neighbors 
to site («, j) and a^t^t' is the lattice vector between sites k and 
kf. Notice that this local vector current does not necessarily 
fall along any of the bonds coming out of the site (nj). 



6.3 Recursion for non-equilibrium Green's functions 



5. Once the left-to-right and right-to-left Green's functions 
are know, the exact inter-slice lesser Green's functions 
are obtained using the Dyson-Langreth equations lISTI 



< , r-L,< 



^n-\ji — ^n-l,«-l U„-l,„G„„ + G^lj „_j t/„_l,„G„„ 



(122) 



In order to evaluate the local current distribution (II lit in the 
most general case, one needs to determine the exact lesser 
Green's function, G^. The procedure is the following. 

1 . We start with the retarded Green's functions of the leads 
(which are assumed to be in equilibrium) and use the 
fluctuation-dissipation relations 



gt^-ifLisl-gl), 



and 



8R^-ifRi8R~8"R), 



(113) 



(114) 



where the advanced Green's functions obey g^l ~ (§£) 
and g'ji = (gg) . Here, fi and /« denote the Fermi distri- 
butions in the left and right leads, respectively. 
2. For the left-to-right sweep, we first determine the re- 
tarded Green's function at the «th slice and then obtain 
the lesser Green's function using the expression 



y^L,< /^^i'' yL.< f^L^a 

^n.n ^n.n n.n ^n.n^ 



(115) 



where the self energy due to the coupling of the «th slice 
to all other slices to the left is given by 



^n.n — Un.n- \G , , Un- \ „ . 



(116) 



^n,n+l ~ ^n,n^n,n+l^n+l,n+l + ^«,n ^n,n+l ^n+l.n 



+ 1' 



(123) 



^n,n-l — G^n.n f^n,«- 1 G^„-l.n-l +G„,„C/„.„-l G„^i^„_i, 



(124) 



and 



< , nR,< 



^n+l,n ^ ^n+l.n+l'^n+l,n^n,n^Grn+l,n+l'^n+l,nG„„ 



(125) 



6.4 Current conservation 

Notice that the current through consecutive slices is con- 
served. Let us prove that using Eq. lllll and writing 



(n-l)^n 



- / dETr (G<„_i f/„_i,„ - f/„,„-i G<_i,„) 



(126) 



Notice that G„,'," = ( G„;,^ j andt/„,„_i = ([/«_!,„) . Thus, 
we obtain the recurrence relation 

Gn,n = {G„lnU„,n-l) G,,!^, „_i (G,,;^, f/„,„_i) . (117) 

3. For the right-to-left sweep, we apply instead the analo- 
gous expression 



fR,< _ f~'R,r yR,< f~'R,a 



where 



^n,n — Un,n+l G^^^j ,,^[ t/„+l,„. 



(118) 



(119) 



Thus, we obtain the other recurrence relation 

G^«,';f = (G^;i,';if^«,«+l) G^„+l,„+l (G^;i.nf^n,"+l) • (120) 

4. In order to join the two sweeps to obtain the exact lesser 
Green's function at a given slice, we simply combine left 
and right self energies and evaluate 

Gn.n = G'^i.n \^i)f + ^n.',f) G",„ (121) 

using Eqs. ( 11161 ) and ( 11191 ) for the self energies. 



and 



I,i^(,i-\-i] — —— / '^£Tr(G„, 1 „t/„„+i —U„^]_„G 



';!^(;i+l) 



n+l,n^n,n+l '^n+l,n^n,n+l J 



(127) 



where the traces indicate a sum over all sites in the slices « 
and « + 1, respectively. Since 



r yL.< I /-•< yL.a 






J J (~'< yL,r ^< I yL,< f^a 



r-^< JJ f^r yR,< I y~^< yR.a 

and 

jT /~^< yR.r f~^< , yR,< ^a 



(128) 



(129) 



(130) 



(131) 



it is straightforward to show that Ii„^i)^„ = 4^(n+i)' which 
guarantees current conservation. It is also possible to show 
that the total current leaving any site is zero. 
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7 Dephasing 

The model Hamiltonian (JUi describes electrons within the 
single-particle approximation, disregarding the effects of ele- 
ctron-phonon and electron-electron interactions. Let us call 
as "environment" all the degrees of freedom that couple to 
electrons in the real physical device. Even at very low tem- 
peratures, due to the interaction with the environment, the 
quantum interference between different electronic paths typ- 
ically fades away at lengths larger than the scale £^. The lat- 
ter is called coherence length or dephasing length. Currently, 
£(p in graphene experiments can be as high as few microns at 
low temperatures, decreasing with increasing temperature. 

This Section shows how to incorporate dephasing into 
the RGF method. We follow the phenomenological approach 
pioneered by Biittiker ll50l and D'Amato and Pastawski ifSTl 
and introduce dephasing in the calculations by adding a set 
of voltage probes to the system. These voltage probes act on 
selected system sites and have their individual chemical po- 
tentials adjusted as not drain or inject any net current. The 
voltage probes give rise to dephasing because the drained 
electrons are not phase coherent with the ones injected back. 
The dephasing length £(p is related to the number of voltage 
probes and the strength of their coupling to the system. We 
focus on the linear transport regime and assume that all scat- 
tering within the voltage probes, albeit incoherent, is elas- 
tic. Therefore we neglect any "vertical flow", as defined by 
DattaEl. 

The basic linear response equations are 



1=1 

Ir == .^RR IJ-R - .^RL Ml - X] .^Ri lii , 
(=1 



where Hr^l) denotes the chemical potential in the right (left) 
lead and /i, is the chemical potential of the /th voltage probe, 
i= l,...N,p, where N(p is the total number of voltage probes. 
The Unear transport coefficients ^ need to be determined 
(see below). The left, right, and probe currents II, Ir, and /,, 
respectively, are defined as positive when they flow into the 
system. Since there is no net probe current, we set /, = for 
all probes. Therefore, 



In addition, notice that when all the chemical potentials are 
identical, all currents should vanish. Thus, 



1=1 

Y,^Ri + ^RL-^RR=0, 

'■'=i('VO 



Using Eqs. (11351 ) and ( 1138b . we can then write 



(136) 
(137) 
(138) 

(139) 



£ W.'^ ill,, - PLr) = ^iL (Pl - ^Ir) , 

where Wu ^ =S? and #<,•/ = -B^ if / 7^ /'. Whenever the 
matrix W is invertible, we obtain 



Substituting Eqs. ( I137l i and (1140b into Eq. ( 1133b . we have 

/ff = J«L(jUff-ML), (141) 

where 



^RL = ^RL^ £ ^Ri 



'-1^ 



,■;/ ~^i' 



i'L- 



(142) 



(132) 


\e\^RL-- 


^'^RL, 




\e\^R.-- 


= ^fli, 


(133) 


\e\^iL-- 


= %L, 




\e\^r- 


= ^r, i^i' 


(134) 


\e\^r- 


'■'=1 ('¥'■) 



^l M/ = E ^l> ^i' + ■^iL liL + ^iR Mfl • 



(135) 



Notice that we have expressed the right-lead current (which 
is equal to minus the left-lead current) in terms of the differ- 
ence between the chemical potential in the leads. 

The linear transport coefficients entering in Eq. (1142) 
can be obtained from the conductance matrix of the system: 

(143) 

(144) 

(145) 

(146) 

(147) 

(This makes ^-Jj #^,/ 7^ 0, thus y^ is in principle invertible.) 
The coefficients /^ are positive. This is consistent with the 
assumption that currents run from higher to lower chemical 
potential. The conductances are calculated using the Caroli 
formulas 

^RL = ^Tr [rRC-^+i.o^LG^m+i] , (148) 

%. = ^Tr [r«G;v+i,„,^G;;,,^+i] , (149) 

^x = ^Tr[l]G,';,,oJlGg,,,], (150) 

^^"" [r,G:_^„^,r,G:_,J. (151) 



h 



-Tr 
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Here, we use the pair {riiji) to denote the site coordinates 
of the ith voltage probe. The level width probe matrix Fj can 
be obtained much in the same way as Fr and 11, namely, 
from the probe surface Green's function (see Sec. |2). Notice 
that in the absence of magnetic fields, all cross-conductance 
matrices are symmetric. In addition, since the coupling ma- 
trices Fr, Fl, and Fj are all Hermitian and positive, one can 
easily show that these conductances are all real and positive. 
In practice, we will assume that the set of lattice points 
attached to voltage probes is sparse, so that N,p ^ NM and 
the computation cost of evaluating the cross conductances 
is not too high. The number of voltage probes N(p and the 
magnitude of the Ffs (see Sec. 17.2b determine £,p. 



leading to 



\E-lii\>2t, 



'^^'^-[M',^F^W^.iE-,,<2,. <'^^> 



The retarded Green's function of a slice which is iso- 
lated from other slices should include the self energy of any 
attached voltage probe: 



?n']Uj') 



E-h„iJ,j)-E[iE)J = f, 



,E-h„{JJ') + iQ+, j^f, 
when n, = n and j = ji. 

13 Recursion for dephasing Green's functions 



(155) 



7.1 Self consistency 



The Green's function entering in Eqs. (1148b to Eq. (1151b 
are the exact ones. Therefore, they have to take into ac- 
count the coupling to the voltage probes and, consequently, 
should depend on the chemical potentials {/i,}. These, how- 
ever, depend on the left and right chemical potentials and 
on the cross-conductance matrices £2, which are given by 
Eqs. (1144b to (1146b . Thus, one can see that this calculation 
needs to be implemented self-consistently. There is one triv- 
ial case though, namely, when the bias across the system is 
zero and ^Lr = jXi (equilibrium condition). In this case, Eq. 
(1140b shows that all /x, = jXr and the self-consistency can 
be trivially satisfied. Nevertheless, one still needs to include 
the self-energy of the voltage probes in the local Green's 
function of each slice during the recursive calculation; h„ — > 
h„ +Ei, when / e n. 



7.2 Voltage probe self energy 

The voltage probe partial width I] is related to the voltage 
probe self energy Ej in the standard way: I] = i{^i ~ -Sf ). 
When each voltage probe is attached a single site, f; and E^ 
are just complex numbers and the insertion of the probe self 
energy into the calculation is substantially simplified. For 
a one-dimensional semi-infinite chain with hopping matrix 
element f coupled to a system site through a hopping am- 
plitude tip, the retarded Green's function reads (see Sec. 14.1) 



g'.iE) 



E-lij 

2f2 



1 



4f2 



{E-m? 



(152) 



Thus, the retarded surface self energy of the /th probe is 
equal to 



Em-tU':{E)=['^) ^ 



1-Wi- 



4f2 



{E~lld\ 
(153) 



In order to evaluate the cross-conductance matrices that en- 
ter in Eq. ( 1143b . it is necessary to evaluate the exact Green's 
function between two slices carrying voltage probes, say, ni 
and «2; for instance, see Eq. (1151b . This is a computation- 
ally intensive task and can only be carried out in a rela- 
tively efficient way if all left-to-right or right-to-left local 
Green's functions are stored during sweeps and the exact 
local Green's functions at either n\ or «2 have already been 
calculated and stored. 

Suppose that we want to evaluate the retarded Green's 
function G,,, „,, with n\ < n2- Here are the steps; 

1. Starting with G^^ „j and G^ ^j „ _|_j obtained during the 
left-to-right sweep, evaluate 



^ni,ni+l ~ ^ni.ni '-'ni.ni+l G„j + i „| + i. 



(156) 



2. Next, use the resulting Green's function and the stored 

<5«,+2,,„+2 to evaluate 

^ni,Hl+2 ~ G^ni,Hi+l ^ni + l,«l+2G„j_|_2,ni+2- (157) 

3. Repeat this procedure until G^ ^ _j is obtained. 

4. From the previously calculated G„, „,, evaluate 

Likewise, for «i > n2, we follow these steps; 
1. Starting with G^^ „j and G^ _[ „ j obtained during the 



right-to-left sweep, evaluate 



1-1 



— G„, „, U„,_„,-] G 



'ni.nj ^ni.ni 



ni,ni-l ^n,-l,n:-l- 



(159) 



2. Next, use the resulting Green's function and the stored 

G^,-2,«,-22 to evaluate 

^ni."l-2 = ^ni."i-l^ni-l.«l-2G^«i-2.ni-2- (160) 

3. Repeat this procedure until G^ „ +i is obtained. 

4. From the previously calculated G„, „,, evaluate 

^n[,n2 ^7j[,H2 + l '^'^2+l'"2 ^ri2,n2- v^^^J 

Since no inversions are required in these steps, the cal- 
culation scales as (9(|«i — yi2\M). 
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8 Disorder 

Disorder is ubiquitous in graphene samples, even in those 
synthesized with state-of-the-art technologies. Depending on 
the synthesis method, charge density inhomogeneities 1521 
or substrate irregularities [l53l . intrinsic and extrinsic rip- 
ples 15411551 . strain fields 1561 . surface molecular adsorption 
l57l . vacancies ll58l . and irregular edges ll59l . are unavoid- 
able. The effects of these different kinds of disorder on trans- 
port in graphene have been addressed by several reviews |lL| 
I2l l60ll6n without exhausting the subject. 

The RGF method is flexible enough to address any of 
the above-mentioned types of disorder by a suitable choice 
of the hopping f,y and the local potential VJ in the Hamil- 
tonian (JUi. Since it is based on an atomistic basis, the RGF 
method is ideal for studying numerically short-range disor- 
der effects, whose typical range is of the order of the lat- 
tice spacing. Although not optimized for that purpose, the 
method can also be used to address long-range disorder BO] 

[lOllSlSl. 

In this Section we present some common disorder mod- 
els for graphene and discuss their implementation within the 
RGF method. 



only about the relative magnitude of the potential fluctua- 
tions, 5V /t, but also about the scatterers' range and concen- 
tration: A simple calculation yields fi40i 

If we now recall that there are two inequivalent atoms per 
hexagon and Ahex = "\/3flo/2, we find that 

2 / B \ 4 



^0 



6Ak- jKmv f^y 



flo 



(165) 



9^/3 ^ V t 
where the numerical prefactor is approximately 40.5 and ,yV 
is the total number of lattice sites. 

In the continuum Hmit and using Eq. ( I163I I together with 
the Born approximation (BA), one finds that the transport 
mean free path away from the Dirac point is given by 

2 If 



C- — 



(166) 



where Xp is the Fermi wavelength in the graphene sheet ||3] 
(Af < ll^ for the BA to hold). 

8.2 Off-diagonal disorder- strain 



8.1 Diagonal (scalar) disorder 

Charge density and substrate inhomogeneities can be mod- 
eled by adding a local disordered potential C/(r,) to the lat- 
tice sites in the sample region. One of the simplest models 
for t/(r,) is constructed as follows: We take o4fmp random 
lattice sites {Ra} uniformly distributed as centers of Gaus- 
sian scatterers with a random amplitude Uk taken from a uni- 
form distribution over the interval [— 5y, 5V] . This results in 



.-^n, 



k=\ 



(162) 



where B, is the range of the potential Il3ll401 . The concentra- 
tion of scatterers is «inip = ,-Aim-p/ ■^^ where s/ denotes the 
total area of the sample. 

-1/2 

In the limit of a low concentration of scatterers, «; ' ^ 
^ , the magnitude of the disorder fluctuations is characterized 
by the dimensionless parameter Ko, which is defined from 
the impurity potential correlation function. 



{U{r„,)U{r„,j)) 



KoiJiv) 

2;r^2 



■nj -■/.',/ 



V4«- 



(163) 



where v = ^/3aot/2h is the Fermi velocity and (• • • ) stands 
for the average over disorder realizations. It is easy to see 
that {U{r,jj)) = 0. We note that Kq contains information not 



Let us consider the situations where the graphene sheet is 
subjected to strain. We now address the case where the strain 
field modifies the carbon-carbon bond lengths, postponing 
to the next subsection the discussion of bond distortion due 
to curvatures. 

Modifications in the bond lengths lead to a hopping renor- 
malization. In the Slater-Koster scheme, the carbon-carbon 
hopping term can, in principle, be obtained from the depen- 
dence of the Vppj: on the inter-orbital distance. In practice, 
one relies on semi-empirical parameterizations, such as 

,-3.37(V«-l) 



Vppnil) 



-te 



(167) 



where I is the bond length and a is the inter-atomic distance 
in the honeycomb lattice. The decay rate is adjusted to fit the 
experimental result dVpp-n,/dl = —6.4 e'V/A. With the help of 
Eq. ( 1167b . local bond length deformations 8lij = lij — a are 
translated into changes in the hopping integrals. The latter 
are easily accounted for by the RGF method. 

The macroscopic theory of elasticity can help to trans- 
late the strain field acting on a graphene sheet into modifica- 
tions of the hopping integrals. That is because the tensions 
along the graphene membrane change very slowly on the mi- 
croscopic scale. Hence, the changes in the bond lengths 5lij 
can be approximated by a smooth function 5l{x,y), where 
{x,y) is the position of the ith site in the honeycomb lattice. 
In turn, the 5l{x,y) can be related to the strain fields by the 
elastic theory Il65l . Reference Il64l . for instance, writes the 
strain tensor e for the case of uniaxial strain and shows how 
to relate e to the bond length deformations. 
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8.3 Off-diagonal disorder - ripples 

Let us assume that there is a ripple structure in the graphene 
sheet. The ripples can be described by a scalar field h{x,y) 
which represents the out-of-plane displacement of the car- 
bon atoms at a given location {x,y). A non-homogeneous 
h{x,y) ripple-field modifies the atomic orbital overlaps and, 
hence, the hopping terms in the tight-binding model. 

Neglecting bond length stretching due to strain, the rip- 
ples affect the nearest neighbor and next-to-nearest neighbor 
hopping matrix elements: r,y = ?■ • + 5f,-.y, where t^- is the 
hopping between sites / and J in the absence of ripples and 



5t,j 



-E,j[{uij-V)Vh] 



(168) 



M 



and 



V, 



ppa;ij 



-2.7eV, 
-0.1 eV, 



5.8eV, 
1.4 eV, 



for n.n. 
for n.n.n., 



for n.n. 
for n.n.n. 



(169) 



(170) 



The free energy associated to the ripple field h{r) for 
graphene on a substrate is given by the Gaussian (elastic) 
form I 



(171) 



F = ]- I d^r [K[V^hf + Y[Vhf + v{h-sf] , 



where K is the bending rigidity, 7 is the interfacial stiffness, 
and V is a coupling constant for the pinning of the graphene 
sheet by the background roughness s{x,y). Typically, JC w 1 
eV. The other two parameters, 7, and v, will depend on the 
substrate. 

In order to generate the appropriate random /j(r), we will 
assume that the temperature is low and thermal fluctuations 
can be neglected (this hypothesis could in principle be re- 
laxed). In this case, the fluctuations in h{r) follow those of 



the substrate: Minimizing the free energy in Eq. (I17U . we 
obtain 



[kV'^ - 7V^ + v) li{x,y) = vs{x,y). 



(172) 



Using the Fourier decompositions li{r) = Jlqh{(l)e"*'^ and 
s{r) = Lq'f(Q)^"' "" we obtain 



(K-^'* + 7^" + v)/i(q) = vi(q). 
As a result, 

v^ {S{qi)s{q2)) 



{h{qi)h{q2)) 



( «■?! + m + v){ Kql + 7^2 + ^) 



(173) 



(174) 



If the background has white-noise fluctuations, (i(qi ) ^(q2) 
i^5P)(qi+q2),weget 



(Ji{qi)h{q2)) = 



{Kq\ 



7?i + v)2 



Since 



(/.(r)/z(0))= £ {h{q,)h{q2))e'^^-\ 
we finally arrive at 



{h{v)h{Q)) =1^- 



,.2 J. 



„'qr 



(175) 



(176) 



(177) 



Here, u,; is the unit vector connecting sites / and j and Eij — heir) = V 

hi /3 + yppa;ij/'2; where Vppajj describes the overlap of the •• 

(7-orbitals (the effect of the a orbitals is only negligible in 
the absence of bending). The scales are the following: 



We now need to find a way to generate membrane profiles 
li{x,y) that satisfy this correlation function. The solution is 
simple: Let us introduce 



„'(q-r^ 



^T + 7T + ^' 



(178) 



where the phases {^q} are uniformly distributed in the inter- 
val [0:27r) and uncorrected, except that 0q = — ^-q in order 
to define a real he. It then follows that 



2 2 



^"'■^'^"^■^'^^ "qS. {Kq\ + Yq\ + v){K4 + yql + v) 



'■qi-i- /g'(<l'qr 



^i^q"' 



,.2,2 



-7q'2 + v)2 



„'qr 



(179) 



which is exactly equal to Eq. ( 11771 ). 

Equation ( 1178b needs to be adapted to a strip geometry. 
First, it is clear that, in rectangular coordinates, 



hc{x,y) = Y. L 



IVSQ 



qx qy>Q 



Kq^ + yq 



2. 



(x^.v+y^v + ^^^g,,) 



(180) 



Moreover, we can assume ^^ = Inux/Nx and qy = iKUy/Ny, 
with n., = -A?,/2, . . . , (A^, - 1 )/2 and «,. = -Ny/2, ...,{Ny- 
l)/2, where A'x x A^y is the number of grid points in real 
space. One could write the real space grid using the primi- 
tive lattice vectors of the hexagonal (actually triangular) un- 
derlying system. However, since the strip has a rectangu- 
lar geometry and the field h{x^y) is defined over a coarse 
grained lattice (hydrodynamic continuum limit) which does 
not need to reflect the underlying atomic structure. Thus, a 
rectangular mesh suffices and we can rewrite Eq. (1179) as 

Nx-l 'Vy-1 

2\'SQ 



hc{x,y) 



y y 

A'v„^o'^^^ + 7'72 + v 

«v= — y- y 



2k X hy^ 



- ^«v.nv 



(181) 



with q^ = {2n)'^[{njc/Nxf + {ny/Ny)\ where we have im- 
plicitly assumed Ny to be odd. 
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8.4 Vector potential (off-diagonal) disorder 

The off-diagonal disorder discussed in the previous Sections 
can be cast in terms of a random vector potential, as nicely 
reviewed in Ref. 1671 . Here we briefly present the mapping 
of the tight-binding hopping disorder 8tij into a random vec- 
tor potential A{x,y). We then show how to implement this 
kind of disorder in the tight-binding model. 

The continuum limit of the honeycomb lattice tight-bin- 
ding model near the neutrality points translates into a Dirac 
equation. Using the Bloch states of one (A or B) of the tri- 
angular sublattices (i.e., the underlying Bravais lattice) that 
constitute honeycomb lattice, one obtains the Hamiltonian 



Hab 



1+e 







1 



„-(k 39 



.ikai 







„!ka2 



(182) 



where ai = (flo/2,flo\/3/2) and 32 = (— flo/2,flo\/3/2) are 
the primitive (Bravais) lattice vectors. Here we use the lat- 
tice vector conventions of Ref. Uj. The Hamiltonian H acts 
on a spinor whose components are the envelop wave func- 
tion amplitudes at the sublattices A and B. The correspond- 
ing low-energy dispersion relation shows two inequivalent 
cones (valleys) centered at the ^-space points 



K 



4-K 
3flo' 



,0 



and 



3fl() 



(183) 



Expanding k around these points, one obtains the Hamilto- 
nians 






k, - iky 
k, + iky 

-k^ - iky 
-k, + iky 



vn{k,d_, + kydy) (184) 

= vh{—k^a^ + kyay), 

(185) 



where kx and ky are measured from the cone vertices. The 
real-space, continuum version of the Hamiltonians can be 
obtained by replacing kx with —idx and ky with —idy. As a 
result. 



Hk = vh 
Hk' = vh 



-idx - dy 
'idx + dy 

idx - dy 



idx + dy 







(186) 



(187) 



When a vector potential is present, the substitution is instead 
kx^y — > —idx.y + ^Ax,y (with e > 0) and the Hamiltonians 
become 



Hk^vTi 



Hk' —vh 







-idx-dy + -§-^{Ax-iAy 



-idx + dy + f(Ax + iAy) 







(188) 



dy-MAx-iAy) 



idx-dy-§XAx + iAy 



tic' 





(189) 



We now show how a local (long-ranged) distortion in the 
lattice gives raise to Hk(k') ^s above. Let us assume that the 
three nearest-neighbor hopping amplitudes of any given site 
are not equal: Calling them Iq, t\, and t2, we have to rewrite 
Eq. ( I182l i in the form 



^^^ ' tQ + tie-'^-^^+t2e-'^-^^ 



fo + fie*'^'+f2e*''- 



(190) 



Expanding around the same K and K' points, and assuming 

|fi-f2|< |fi+f2|,wefind 



Hk = vh 



kx - iky 



kx + iky 
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Hk' ~ vh 
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Thus, we can define two vector potentials, one for each cone 

h+t2 



^^-VeV^- 2 



A^ = -^(fi-f2), (193) 
ve 2 



iK 



and 



fo- 



fl+f2 



c V3 
ve 2 



{ti-t2). (194) 



Notice that A^ = — A^ , as expected from time-reversal sym- 
metry considerations. 

We can use the vector potential as a gauge field that pa- 
rameterizes local fluctuations in the nearest-neighbor hop- 
ping matrix elements. For this purpose, it is useful make a 



single-cone approximation and rewrite Eq. (|1 93) in the form 


2ve 
to{nJ)=t + —Ax{nJ), (195) 
c 


/ .X ve 


AxinJ)+ Ay{n,j) 


(196) 


/ .X ve 


AxinJ) A,,(«,;) 


(197) 



where {n,j) are the coordinates of sites belonging to one of 
the sublattices. Notice that this type of off-diagonal disorder 
does break particle-hole symmetry. 

Another way to proceed and get the same results is to use 
a Peierls substitution in the hopping matrix elements 
such that 
e'^-ai . e'(k+j^A)-ai 



„'ka2 



,<-(k+j^A).a2 



(198) 
(199) 
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in Eq. (1182b . For instance, expanding on A and around the 
K point, we obtain 



1 +e'(''+7fe'^)-''l +e'(''+j|:A)-a2 



hv [k^ - ik^) H (Av - /Av) . 

c 



Then, setting 



to 



ti+ti 



~iV3- 



— (A., - iAy) , 

c 



(200) 



(201) 



2 2 

we arrive at Eqs. (|195t to (|197| i. 

This concludes the demonstration that a hopping distor- 
tion can be mapped into a corresponding vector field. From 
the point of view of a tight-binding modeling and the RGF 
method, it seems simpler to model strain with renormalized 
hoppings, thus avoiding issues related to projections onto 
the K and K' cones to preserve time-reversal symmetry. 

The random gauge potential model is also interesting 
for other reasons. It has a physical realization in rippled 
graphene subjected to a strong parallel magnetic field ||69l 
and has been analytically addressed by several authors, e.g. 
Ref. E). 

Let us model the random vector potential by assuming 
that the gauge field has Gaussian fluctuations, such that 



{Aa{n,j)Ap{n'j'))^X8, 



'ap 



„-|r«,;-r„/,/IV2?- 



(202) 



where X measures the strength of the fluctuations and ^ is 
their correlation length. One way to generate this correla- 
tion function is to define at the nodes of a regular lattice of 
constant Og two sets of uniformly distributed random num- 
ber^ {c"}k=i,...,.v and to define the gauge field through the 
expression 



(203) 



Aa{nJ) = 


f -^ 


-R.lV«^ 


where 






k=i 


-|r,,j-R,|V^2 





(204) 



and A = f^lug/^)^ /2% when JV -^ °° maintaining ^Og < 
oo. This construction implicitly assumes that a^ <C ^ . 

There is a useful way to quantify the fluctuations of the 
vector potential. Let us denote (^^) the rms value of the 
magnetic flux piercing a region of area £/ <^^. Then, ( 0' ) r 
£/^ {B^), where B = (9,A,. — dyAx- Following the steps shown 
in AppendixO we find that {b^) = f{aglE,f /2k^^. 

We can now redefine the vector potential to absorb the 
prefactor e/hc: A = {e/hc)A. Likewise, in order to get rid 
of the prefactor in the expressions used to generate the vec- 
tor potential, we introduce / = {hc/e)f and X = {tic/ef'X. 

■' One set of random numbers for each a -component of the gauge 
field, such that (c") = and {cfc"') = Skk^Saa'- 



Then, we can write the following expression for the esti- 
mated rms value of the random magnetic flux in units of the 
flux quantum (^o = hc/e): 



2k- 00 ^ 27r V2^ U' ' "' 



(205) 



which implies X = {5(p)^^^/£/^. Thus, the relation between 
the rms flux phase piercing an elementary hexagon £4ex = 
\/3 flo/2 and the vector potential intensity is 



V 3 \aoJ Gg ' 
or, equivalently, 

/hex _ [j^f^y^ 



(206) 



(207) 



On the other hand, if we set ^ppie = ^ ^ to denote the typical 
area of a ripple, we obtain 

/.„,.. V2J^, |i.^/|faV|5, ,208) 



TT \ ^ y Qg 



8.5 Edge disorder 



Etching a graphene sheet to produce nanoribbons always 
leaves behind some roughness at the edges. When the ir- 
regular shape of the boundaries of the propagating region is 
very pronounced, it leads to the formation of "bottlenecks" 
and "cavities", which tend to increase charging effects and 
lead to Coulomb blockade oscillations of the conductance 
llTOl . However, even mild amounts of edge disorder can af- 
fect dramatically electronic transport in nanoribbons. In this 
case, it has been proposed that for long enough nanoribbons, 
Anderson localization (and thus an insulating behavior) can 
develop 17 1 1172117311 . Insulating behavior, albeit of a differ- 
ent nature, is also expected in "perfect" nanoribbons due to 
lattice symmetric breaking caused by the deformation of the 
chemical bonds involving carbon atoms at the edges, as re- 
vealed by DFT calculations B3]| 

Edge disorder can be simulated by considering slices 
with random numbers sites (see Fig. fl4l i: For instance, we 
can draw the number of sites M„ of the nth slice randomly 
according to the Gaussian distribution 



P{Mn) = 



1 



2%5M 



^-(Mn-uf ll&M^ 



(209) 



where M is the average number of transverse unit cells in the 
nanoribbon and 5M is its standard deviation. Other distribu- 
tions can be investigated straightforwardly. The slices are 
concatenated such that hopping matrix elements connecting 
sites which fall into empty spaces are set to zero (although 
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this can be avoided when programming the recursive calcu- 
lation by using nested loops with variable ranges). Thus, one 
may think of this procedure as a random removal of sites at 
the edges of the nanoribbon. This approach has been used to 
study the existence of localized states in graphene systems 
CD. 

Another approach, which tries to mimic the effect of 
etching, is explained in ff3\. Again, the numbers {M„} are 
considered random variables, but their generation follows a 
different procedure. One visits sequentially each edge site 
(at the top and bottom) and elects to removes it (or not) ac- 
cording to a probability pi. Certain sites, when removed, 
require the removal of neighboring sites as well, as an edge 
configuration where carbon atoms have a single bond are 
not stable (unless both dangling bonds in the carbon atom 
are pacified, but this is not likely to occur during etching). 
After this first sweep of edge sites, a second sweep follows, 
but now sites are removed with a probability p2- One can 
continue repeating this procedure, using a different removal 
probability at each sweep, until the desired amount of rough- 
ness is obtained. 



9 Some Numerical Results 




Fig. 10 Results of linear transport calculations for clean graphene 
sheets: Conductivity (defined as ff = LG/W) and Fano factor as a func- 
tion of the Fenni energy in the contacts {Vg = 0). The band in the square 
lattice leads was offset such that its middle coincides with the neutrality 
point in the graphene sheet in order to increase the density of states at 
the contact and thus mimic a metallic lead. Armchair edges, M = 360 
and W = 70 (aspect ratio W/L = 5.2). 



In this Section we show some representative results obtained 
with the RGF method. We begin by addressing the case of 
ballistic transport in graphene sheets where analytical results 
are known and serve to benchmark the numerical method. 
Next, we discuss the case of graphene sheets with long- 
range disorder, where the RGF method was used to clarify 
the controversial issue of the "universal conductivity mini- 
mum". Finally, we present results for graphene nanoribbons , 
showing how the method can be used to calculate LDOS and 
the local current density. 



9.1 Ballistic transport in clean samples 

Here we present results obtained for the case of ballistic 
transport in graphene sheets. We consider mainly the arm- 
chair orientation, since this, in the clean limit, provides a 
band structure and dispersion relation very similar to a quasi- 
one-dimensional projection of the Dirac fermion model and 
is more suitable for scaling analyzes. 

First, in Fig.jTO] we show results for the clean limit (no 
bulk or edge disorder) for a short ribbon, keeping the back 
gate voltage fixed to zero (neutrality point) and varying the 
Fermi energy in the contacts (zero bias). The numerical data 
is compared to the analytical expressions derived by in Ref. 
||39l . The agreement is quite good for large systems and be- 
comes worse when the system is too small (not shown). In 
particular, a strong asymmetry and a lack of well-defined 



oscillations occurs if the system is not large enough (not 
shown). 

The ability of the recursive method to get precise results 
for clean systems is clear also in Fig. [TT] where the calcu- 
lations are performed for different aspect ratios. The devia- 
tions from the analytical curve only occur when the system 
is too short and evanescent modes dominate transport. 

In Figs. [12] we show the conductance and the current 
density in the linear regime for small flakes with armchair 
and zigzag edges. Here, the leads are also honeycomb lat- 
tices. In this case, the conductance steps can be easily un- 
derstood from the energy dispersion relation of graphene in- 
finite ribbons 11741 . namely, the dimensionless conductance 
G/{2e-/h) is given by the number of bands crossing the 
Fermi energy E. 

There is no such simple explanation for the current den- 
sity. Notice that the notable difference in the current distri- 
bution for armchair and zigzag orientations at £ = 0: While 
in the latter the current is primarily carried by edge states, 
in the former the current is uniformly distributed across the 
flake. As one moves a just little bit away from E = 0, the 
current distribution for the zigzag flake changes drastically, 
with nearly no current running at the edges. This result is 
related to the the fact that for zigzag nanoribbons the E =0 
states are strongly localized at the edges. As soon as l^l > 
both edge states and edge currents disappear, even in the 
case of a single conducting channel. For the armchair orien- 
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Fig. 11 Results of numerical simulations of clean ribbons: conduc- 
tivity and Fano factor as a function of the Fermi energy in the con- 
tacts (Vg = 0) for different aspect ratios, (a), (b), (e), (f): Square lat- 
tice contacts; (c) and (d): Armchair honeycomb contacts. For all plots, 
M = 120. In plots (a), (b), (c), and (d), the value of W are 12 (dashed 
line), 24, 36, 72, 96, 108, 120, 148, and 200 (dashed-dotted line). The 
thick solid line corresponds to the analytical result 1391 . 
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Fig. 12 CuiTent densities (a,b,d,e) and lineai' conductance steps (c,d) 
of two small ballistic graphene ribbons. In (a,b,c,d) the arrows repre- 
sent the current densities (in arbitrary units) evaluated at different sites 
using Eq. fTTH . Armchair edges, M = 12, W = 20: (a), (b), and (c); 
zigzag edges, M = 12, A' = 21: (d), (e), and (f). Energies: £ = for 
plots (a) and (d); E = 0.3f for plot (b); E = O.Olf for plot (e). ampli- 
tude. 



tation, the change in the current distribution for increasing 
energies is less drastic. 

It should be stressed the edge current densities of zigzag 
nanoribbons change both quantitative and qualitative if one 
switches from nearest-neighbor ll75l to next-nearest-neighbor 
tight-binding models 1761 . The issue of which model is ap- 
propriate is tied to the desire to fit DFT calculations B31 or 
to explain experimental manipulation and characterization 
of nanoribbon edges llTTl . 



9.2 Disordered graphene sheets 

The conductivity minimum Ob observed at the charge neu- 
trality point in graphene monolayers has been a subject of 
intense debate, which is reviewed, for instance, in Ref. ||2l . 
Here, we show how the long-range Gaussian correlated po- 
tential can be used to investigate the value of CJo-u 

** The discussion and results that follow complement the material 
presented in Ref. J2]- 



In the diffusive regime, in general, the system geometry 
has little influence on the transport properties which allows 
one to express the average conductivity as a = (L/W)(G), 
where L is the system length and W its width. We use the 
same setting as in the previous subsection, including now 
long-range Gaussian disorder in the device region. To gen- 
erate the data shown in Fig. [13] four different aspect ratios 
were considered as well as several values of Kq and ^ja. 
The average conductivity Oq obtained from (G(Vg = 0)) is 
plotted versus L scaled by t . The parameter i* depends on 
/iTo and ^ . We identify i* with the elastic disorder mean free 
path I. 

Let us summarize the results shown in Fig. [13] Two clear 
regimes can be identified. For Lji <C 1 , the probability of an 
electron being scattered by disorder as it traverses the sam- 
ple is very small. This corresponds to the ballistic regime, 
where scattering occurs mainly at the sample edges and trans- 
port properties are dominated by the sample geometry. Note 
that when L/£ < 1 , cTq approaches the prediction for the pure 
ballistic case [391 . indicated by the arrows in Fig. [13] In con- 
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Fig. 13 Conductivity minimum 0b in unit of 2e^//i as a function of 
the system size L scaled by the electron mean free path I. The results 
correspond to the average over 10^ • • • 10* disorder realizations for L 
ranging between 50 and 500 Qq- The colors represent different aspect 
ratios WjL. The symbols stand for the values of the dimensionless dis- 
order strength A'o . The arrows indicate the analytical value of the con- 
ductivity minimum in the ballistic limit |39| , which depends on WjL. 
The dotted line gives the diffusive ln(L/£) behavior (2. 



trast, when LjiL ^ 1, the system becomes diffusive and ge- 
ometry affects transport weakly. Figure [TSfclearlv shows this 
crossover For the diffusive regime, Ljl ^ 1, the conductiv- 
ity is proportional to ]n{L/£), in agreement with the non- 
linear sigma model prediction H. The mismatch between 
the numerical prefactor for the logarithm and the value char- 
acteristic of the symplectic class may be related to the finite 
contact resistance [621 present in our simulations. 

These simulations suggest an explanation for results ob- 
tained in transport experiments at the charge neutrality point. 
In the coherent diffusive regime, the conductivity minimum 
has significant sample-to-sample fluctuations and its average 
shows a weak (logarithmic) dependence on the mean free 
path. Typical diffusive experimental samples have L/i w 
1 — 10 and (7o ~ A-e" /h, similarly to what is shown in Fig. [13] 



9.3 Nanoribbons 

For nanoribbons, both bulk and edge disorder play a role 
in electronic transport. In the absence of band gaps, long- 
range disorder does not suppress conductance significantly 
and a perfect conducting channel exists near the neutrality 
point Il78ll79|[80l . The story is quite different for short-range 
disorder Bulk imperfections (lattice defects, impurities, or 
adsorbates) and edge imperfections can lead to strong local- 
ization due to backscattering and enhanced destructive in- 
terference 07211731 . To illustrate this point, Fig.[T4]shows the 
rapid smearing of the linear conductance steps of a nanorib- 
bon when even a small amount of edge sites are randomly 
removed (i.e., etched out). This shows how challenging it is 
to observe conductance quantization experimentally in these 
systems. 
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Fig. 14 (a) Linear conductance of edge disordered nanoribbons with 
armchair edge orientation, M = 18, A' = 200, averaged over 100 real- 
izations, as a function of energy. Only one etching sweep is used, but 
results for three different values of the site removal probability p\ are 
shown, (b) Typical realizations used in (a) for value of pi considered. 



The appearance of localized states in nanoribbons with 
edge disorder is demonstrated in Fig. [15] where the linear 
conductance and the local density of states for a nanoribbon 
are shown in the cases of perfect and irregular edges. 
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A Steps in the linear conductance 

Let us show that the surface Green's function in Eq. | |67) leads to the 
expect steps in the linear conductance. For this purpose, let us begin 
by noticing that, in the case of a square lattice lead, only propagating 
modes yield a finite level width: For |£ — ev| < 2?;,, 



fv = -2Im [Zv] = 2Im (gv) = 2fv sin 



(210) 



where sin ^y = y'l — (E — EyY /At^, in which case we can write gv = 

In order to obtain the retarded Green's function across the system, 
we add one slice between the left and right contacts and use the fol- 
lowing expression, easily derivable from Eqs. ( I32K ( 1331 , ( I4U and l l42t : 

Go.2 = f;gife'-f?gt)"' (211) 

Since Gqj depends solely on gi, we can rewrite in the propagation 
mode basis, in which case the Landauer formula is reduced to [see Eq. 

^ = ^'fv(Go.2),fv(Go,2);, 



(212) 



24 



Caio H. Lewenkopf, Eduardo R. Mucciolo 




Putting all together, we find that 

^=i:'n^i(Go.2),f 

V 

=^ 1 = # propagating modes for a given E, 

V 

which is the expected result for a clean ballistic system. 
B Peierls hopping phases 
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Fig. 15 (a) The linear conductance of a short nanoribbon (zigzag 
edges, M = 24, N = 91 sites) as a function of energy. The solid line 
corresponds to perfect edges while the dashed line corresponds to the 
edge disorder realization shown in the inset (three etching sweeps with 
p\ = 0.3, P2 =0.2, p3 =0.1). Surface plot of the local density of states 
(arbitrary units, brickwall lattice representation) for the same nanorib- 
bon at the energy value highlighted in (a): (b) clean case; (c) edge dis- 
ordered case. Notice the appearance of localized states that traverse 
the nanoribbon when edge disorder is present. (The actual LDOS was 
convoluted with a Gaussian profile to smooth out high-frequency os- 
cillations.) 



Here we evaluate the phase of the hopping matrix elements between 
any two arbitrary sites due to the presence of a perpendicular magnetic 
field. We pick the vector potential in the generic Landau gauge A^ = 
(a - i)By and Ay = aBx, with < a < 1. The (directional) Peierls 
phase between two neighboring sites k and k' is given by 1681 
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(215) 



where 9^i^/ is the angle that the segment k-k' makes with the x axis. 

Notice that (pi, ji = — (p^r./ j-. If we sum over all the bond phases 
around the perimeter of a hexagon, we obtain £(p = ^ecBa^/2h = 
2lt{<P/4>o), where ff>o = h/ec (flux quantum), and <P = BAhex, with 
'^hex = x/S ^0/2 being the area of the hexagon. 



C Random Flux Estimate 

Let us estimate the nns value of the random magnetic field produced 
by the random vector potential: 



{B^) = {(dAy-dyA,f) 
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